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Stochastic Modeling: Ideas and Techniques 

Donald P. Gaver 



1 . Introduction 

The primary purpose of this chapter is to summarize the 
contents of lectures on stochastic modeling presented at the 
Universite Libre de Bruxelles (ULB) in the period March-May, 

1981. Much of the material selected for presentation was from 
the standard menu of probabilistic topics typical of a second 
course as given to engineers, operations researchers, 
statisticians, or computer scientists. An attempt was made to 
emphasize a modeling attitude rather than details of mathematical 
rigor, illustrating with problems and techniques that are not 
often prominent in such courses. For example, attention was 
given to problems of, and models for, redundant system reliability 
and availability, queueing with priorities, first-passage times 
and areas under path functions of stochastic processes, (total 
waiting times), and various other topics. Also included was a 
brief account of aspects of modern data analysis, with the 

implication that its usefulness is significant at the pre-modeling 
and model-assessment stage of an investigation. 

A secondary, but gratifying, purpose is to briefly report 
on cooperative work initiated with faculty and students at ULB. 

I wish to mention the enjoyable collaboration with Dr. Guy Latouche 
on development of efficient computational methods for repairman- 
like Markov models in random environments, and with Ph. Collard 
on the application of sculptured distributions in the simulation 



evaluation of certain scheduling algorithms. The interest and 
warm hospitality of Prof. Guy Louchard, head of J the Dept.' of 
Computer Science at ULB, was also much appreciated. 
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2 . 



The Total Modeling Process: Brief Overview 



It is coining to be recognized that the topic of mathematical 
modeling (including stochastic modeling) exists in its own right 
as a subject suitable for a formal university course; see 
Bender [1978], The modeling step is part of a process of several 
stages or steps; these may be expressed as follows (Gaver and 
Thompson [1973] ) : 

(a) Identify the general problem area or situation ; identify 
specific questions concerning that area. 

(b) Obtain and analyze subject-matter information and data 
relating to the problem area. Often an examination 

of such information and data will suggest suitably 
formulated questions, as in (a). 

(c) Construct a preliminary model, or models, representing 
the important features of the situation. Deduce some 
model implications. 

(d) Refer the result of (c) to subject-matter specialists 

and decision-makers for qualitative critique; revise 
the model accordingly. This likely means re-doing 
(a) - (d). 

(e) Assess the empirical validity of the model to the degree 
possible. Check the sensitivity of model conclusions 

to changes in model assumptions (sub-model inputs), 
and to data variations. Submit to judgement by 
subject-matter experts -- but anticipate differences 
of opinion! The modeling, and re-modeling, process 
may help to reconcile such differences. 
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(f) Compute required answers to interesting questions. 

Assess the degree of uncertainty in these answers 
possibly resulting from model mis-specif ication , data 
bias or other deficiency, computational error, and 
sampling error in estimates of basic parameters or 

in simulation results used to supply model implications. 

(g) Communicate, and aid in implementing, the results of 
the model. 

(h) Monitor the situation for possible changes in the 
environment, and hence for the necessity to change 
the model. 

Of course the emphasis in these notes (and in the lectures, 
was) upon the actual modeling step, (c) . However, some attention 
was given to the display of data for pre-model examination 
(Tukey's exploratory data analysis), and to model parameter 
estimation techniques, particularly those robust methods that 
attempt to deal with questions of data deficiencies. 



3 . 



Topics in Outline 



In this section we out-line the basic contents of the 
lectures. These were in general arranged so as to first present 
mathematical definitions and properties, and then to illustrate 
in terms of sample models for various situations. 

( 1 ) Review of Probabilistic Concepts, Particularly Conditioning . 

In this lecture the following basic notions of probability 
were defined or reviewed: random experiment or trial, sample 

or event space, events and combinations of events, probability as 
a function with rules for combination, conditional probability, 
independence, and Bayes' Theorem, random variables and their 
moments or expectations, transforms (characteristic function, 

Laplace transform, and generating function) and their moment- 
generating, convolution, unicity, and continuity properties, plus 
properties of conditional expectations. In addition, certain 
classical univariate distributions were reviewed (Normal/Gaussian, 
log-Normal, exponential and gamma, etc.) 

By way of illustration, a simple problem of equipment (or 
possibly software) unrealiability was considered. 

Situation: Suppose a system is made up of components that 

individually fail after a time because of the action of faults ; 
the latter may be the result of component mis-design, or attributable 
to bad installation or adjustment ("human error"), or to a mistake 
in computer program coding. We wish to relate system failure rate 
to initial fault content. 

Model: N is a random variable (rv) representing the 

number of faults initially installed unwittingly in the system. 

Let (T^, i = If 2, ..., N} be the sequence of rv describing 
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the failure time of each fault, measured from time at which the 
system begins use; Tb may actually be the time at which the 
service of the particular component is first requested. Suppose 
the system fails at 

X N = min { T 1 , T 2 , T^} . (1.1) 

Under very simple conditions, namely that all components have 
the same distribution, F(t), of failure time, and all failure 
times are independent, simple conditional probability arguments 
yield 

p(x H > t) = e n £ (i-f ( t) } N ] (1.2) 

= g N (F(t)) 

where g N is the generating function of the number of faults 
originally sown in the system, and F(t) = 1-F (t) is the 
survival time distribution, per fault. It is easy to see that 
P{X N = 00 } = p{n = 0}, possibly > 0 , so the derived distribution 
of X N is quite possibly dishonest. Note that while in general 
explicit expressions for expectations cannot be obtained (may 
not even exist) , such summaries as the median, 90% point, etc . , 
may if g (f ( t) ) can be explicitly inverted, e.g. for 
N ~ Poisson and F ~ Exponential. The simplistic assumption 
of the model may be relaxed, allowing for different distri- 

butions, dependence, and so on, and an additional random death 
time, D, applying to the total system can be introduced to in- 
duce eventual failure (of physical equipment), or biological 
death in finite time. There will be less analytical tractability , 
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but simulation may be used to assess system behavior. Statistical 
estimation problems may be addressed as well; a suitable version 
of (1.2) will provide a likelihood function. 

Another example of the applicability of a simple conditioning 
argument is the following. 

Situation : When an individual speaks on a telephone or 

telecommunication channel the conversation is an alternating 
sequence of talk-spurts and pauses. Similarly, a job being 
processed on a computer goes through an alternating sequence of 
CPU (compute) times and 10 (input-output) times. Model the 
total time of the conversation or job processing time, and 
particularly the joint distribution of busy and idle segments. 
Model 1 ; Let {X if i = 1,2, ..., K} and {Y^ i = 1,2, ..., K} 

be the durations of talk spurts and pauses, respectively, and 
let K be the number of each. The simplest model assumes 
{X^} and { Y ± } to be independently and identically distributed 
(IID) sequences of rv, and themselves to be conditionally in- 
dependent, given K, also a rv. The joint distribution of total 
talking (or processing) time, X, and total pause time (10 time), 

Y, is thus, by simple conditioning. 



r- k* k* 

P{X < x, Y < y} = l F (x ) • F (y) P{K=k} 

k=l X Y 



(1.3) 



The joint Laplace transform is (s^, s^ , ^ 0) 




oo 



l [F x (s 1 )F y (s 2 ) ] k P { I<= k } 



k=l 



(1.4) 



= g K [F x (s 1 )F y (s 2 )]. 
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where g v is the generating function of K. Put s, = s_ = s 
J\ 1 z 

to recover the transform of X+Y = L, the total conversation 
length. In case 



~ Expon(A) and ~ Expon(y) and K ~ Geom(a), independent 

k 



-s,x -s 0 Y* 


00 


r h. > k 


( \ 


1 2 
r\ o 


= l 

k=l 


A 


y 




[A+sJ 


1p+s 2 J 



(1-a) a 



k-1 



(1.5) 



Ay (1-a) 



( A+s^ ) (y+s 2 ) - Aya 



and 



E[X] = [ A (1-a) ] - 1 , E[Y] = [y (1-a) ] 1 



( 1 . 6 ) 



Furthermore (a = 1-a) 



E[L] = E [x] + E[Y] = (1.7) 

Aya 

and 

Var [ l] = (E[L]) 2 ^ (1.8) 

Aya 

Notice also that the mechanism of randomization of a sum , or 
mixing (see Feller [1966j)which has given (1.3) may be used to 
generate families of bivariate (multivariate) exponential 
distributions for other modeling purposes. 

Model 2: A plausible alternative to the above model assumes 

X^ and Y^ are not independent, being possibly positively 
correlated -- a long talkspurt tending to result in a long 
pause (response by conversationalist) . Most simply, 
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= BX^ , B > 0. Then again transform in the X 
to get 



Expon (A) case 



E 



- S 1 X - s 2 Y 
e e 



A (1-a) 



S 1 + ^ S 2 + ^ l _a ) 



or if B - A/p which preserves the marginals of Model 1, 



A (1-a) 



s. + — s n + A (1-a) 

1 p 2 



Apa 



S l y + S 2^ + 



(1.9) 



It is now immediate that the marginal distribution of 

X ~ Expon(Aa), Y ~ Expon (pa) and that now the df of the 

Apa^ 



total time, L, is simply Expon 



, A+p 



-- a much simpler form 



than that occurring in Model 1 above, which involves a Bessel 
function. The variance of L in Model 2, being (e[l]) , is 

also larger than that for Model 1, see (1.8), suggesting that 
the former model has a longer tail, hence predicting a greater 
proportion of extremely long conversations. 

The above illustrates that the same situation can easily 
give rise to two — or more — different models, depending 
upon the manner in which stochastic assumptions are introduced. 
At best, the introduction should be guided by observed data; 
at least, sensitivity analyses using different assumptions can 
outline the range of specification uncertainty. 
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(2) Models Involving Repeated Trials 



A great many situations may be initially modeled in terms of 
repeated independent trials, where this means that on each of a 
possibly countably infinite number of occasions a trial (or 
experiment, or observation) is performed, with outcome 
(possibly a vector random variable (rv) ) on the ith trial; 

{X^, i = 1, 2, ...} are IID rv. Bernoul li trials are a prime 
example: flip a biased coin indefinitely; let 1^ be one if a 
Head (Success) results, and zero if a Tail (Failure) otherwise, 
and assume that probability of success on any trial is independent 
of all previous outcomes. Equally, X^ may be the winnings on a 
bet at occasion i, with X^ in dollars and either positive or 
negative. Or X^ may even represent the increase or decrease 
in a qommon stock price on the New York Stock exchange, according 
to some observers. 

It is convenient, but somewhat more questionable, to adopt 
the repeated trial model for modeling real operational and 
physical phenomena, yet it is often done uncritically. 

For instance the lifetimes (times between failure) of 
computing equipment are often modeled by IID rv, repair times 
likewise, queueing system inter-arrivals and service times as 
well, inventory demand sizes, sizes of deposits of resources 
(petroleum) as well, ... the list is very long, and observational 
support for these assumptions is usually conspicuously lacking. The 
Attraction of the repeated trials model is mainly its mathematical 
tractability , which leads to elegant and appealing results. 
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Often a brief data analysis in terms of marginal distributions 
of observed seems to provide justification. The simple 

repeated trials model cannot well represent, say, systematic 
daily changes in job inter-arrival times, or numbers of jobs 
per hour, at a computer center, seasonal effects on such computer 
system demand, or the influence of other variables such as the 
introduction of a new class of computer users upon a measure of 
computer loading. Also, the model does not well describe the 
sequence of daily rainfalls in a region, nor many other environ- 
mental variables. Some examples follow in which repeated trial 
models seem initially plausible, but which doubtless can stand 
improvement. 

Situation : A structure is to be designed to withstand (wind- 
generated, or seismic or other) shocks. Question: how long will 

it survive the environment stresses, given its initial strength Z? 
Model : Let the structure have strength Z (suitable units) . If 

the ith year's maximum shock is X^, assume {X^} to be IID with 

df F (x) . Then the time T until structure failure exceeds t 
X 



(t = 1, 2, . 



iff YXj^ < Z, i = 1,2, ..., t, so 



P{T > t | Z} = [F X (Z ) ] t , 



( 2 . 1 ) 



the geometric distribution, 



while if Z itself is regarded as random 



P{ T > t} = E 



Z 



P{T > t | Z} = E Z {[F X (Z) ] fc } 



( 2 . 2 ) 



f (E z ([F x (Z)]}) t ; 
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the unconditional distribution is a convex (probability) combin- 
ation of exponentials. It will resemble an exponential, but often 
has an extended right tail. It is obviously important that the 
condition on equipment stress, Z, be removed at the appropriate 
stage -- at the end. Removal of the condition before each 
"strength test" would be appropriate only if the structure were 
completely repaired or replaced with another having the unchanged 
distribution of Z before each trial. 

Note that assuming yearly maximum environmental events to 
be IID is intuitively plausible, but some statistical evidence 
exists for truly long-range correlation in weather data that 
may call even this assumption into question. 

The Bernoulli counting process is an important special case 
of repeated trials, on each of which either Success (probability 
p) or Failure (probability q = 1-p) occurs; {N (t) , t = 1,2, ...} 

is the number of Successes in { 1 , t ] ) . Times (number of trials) 

t between successive Successes are IID and geometrically distrib- 

k— 1 

uted (P{t = k} = q p) . The number N(t) of Successes in 
fixed time is Binomial (p, t) , and is in turn approxi- 
mately Normal (tp, tpq) as t °°. Of course Bernoulli trials 
describe the outcomes of many other repeated trial situations: 
for instance, the number of jobs submitted to a processing 
facility requiring more than x time units of processing may 
be modeled as a Bernoulli counting process with p = 1 - F (x) , 

F^ being the distribution of job processing time. 

Generalizations of Bernoulli trial situations may be (a) 

to variations of success probability with trial number ("time", 

t 

or "space"): P{X. = 1} = p. . Here E[N (t ) ] = £ p. and 

11 i=l 1 

t _ _ , t 

Var[N (t ) ] = l p.^. 1 tp(l-p), where p = r l Pj_, so variability 
i=l i=l 

is under-represented if trial-to-trial probabilities change, but 
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a Normal approximation to N(t) may still hold. A second generalization 
is to (b) independently randomize p in the Binomial distribution, 
say according to a Beta distribution, creating the Beta-Binomial dis- 
tribution. This has found useful in reliability modeling and 
in Bayesian inference. A third and important elaboration/ (c) , is to 
develop a statistical regression model for success probability p 
based conveniently on the logistic function ; 



a+Bu . 

_e 

p i - a+(3u. ' 

1+e 1 



(2.3) 



here u^ (possibly a vector) represents the influence of other 

factors upon success probability. Given observations of the 

form (I., u.) , where I. = 1 indicates success on trial i, one 
11 l 

can estimate a and 3 (vector) by maximum likelihood; see 
Cox [1969 l]. Generalizations to multiple-category situations are 
possible, and computational methods for parameter estimation 
and model assessment have been devised; see Pregibon [1901] . 

The distribution of the number of counts N(t) in t trials 
of a Bernoulli trials process can be computed by making use of a 
forward equation . Let 



Pj (t) = P{N (t ) = j | N ( 0) = 0}. 

Then 

Pj(t) = Pj (t-1) • (1-p) + Pj_ 1 (t-1) -p;(0 < j < t) ; (2.4) 

on the basis of conditioning on events that have happened up to 
t-1. One can generalize to non-stationary success probabilities 
easily : 
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(2.5) 



Pj(t) = Pj (t-1) • (1-P t ) + P j _ 1 (t-l)p t . 

Initial conditions may be 

P Q (0) = 1, P_. (0) =0, j = 1,2, ... 

A further generalization allows success probability to depend 
upon the number of previous successes; then 

p j(t> = p j(t-i) (i-p j(t ) + - < 2 - 6 > 

and the distribution (t) can easily be computed recursively 
given the success probabilities 

p = P{N(t) = j + 1 1 N (t) = j} . (2.7) 

j 

This is a preview of ideas of Markov chains, to be treated later. 
These expressions are introduced to suggest early that the answers 
to interesting and comparatively complex problems can be directly 
computed numerically (in this case iteratively, starting from 
the initial conditions) . Closed-form expressions such as the 
Binomial distribution are handy, and Normal approximations are 
even handier, but one need not modify the facts merely for the sake 
of convenience. 
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(3) Sums of Repeated Trial (IIP) Random Variables; "Large 
Deviations " 

Models for total demand for physical inventory or for 
facility (computer) time often naturally involve sums of varying 
components, modelled as rv. ; thus total demand from n sources, 
or over n time periods, is 

S =X, +X„+...+X. (3.1) 

n 1 2 n 

If X i is the (dollar) profit in the ith year for some enterprise, 
then a financial measure of success is 

S (r) = l x.r 1 (3.2) 

n • n 1 

1=1 

where r is a discount rate (0 < r < 1) . 



Situation . A computer center experiences varying monthly demands, 
X^ for the ith month. Here are answers to several simple 
questions involving sums of X^s. 

The expected yearly (n-period) demand is 

E[S ] = l E[x . ] , (3-3) 

i=l 



the sum of the expected monthly demands, and also 

n 



Var[S n ] = l Var[x i ] 
i=l 



(3.4) 



provided the X^s are uncorrelated . Importantly, as n 



F g , (x) = P 
n 



S - E[S ] 
n n 

/VarLS J 
L n 



< x 



1 2 

X -~z , 

if 2 dz _ . , . 

% j e = <f> (x) , 

-co /2 t r 



(3.5) 
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i.e. S n becomes approximately Normally distributed no matter 



what the distributions of X^, by the Central Limit Theorem, 
provided the components are all of about the same size 

(certainly if they all come independently from the same parent 
distribution with finite mean and variance) . For smallish n or 
distinctly non-Normal components the approximation is improved 
bv an Edaeworth expansion (Feller [19661), wherein for the ecual 



component example 



Fg ,(x) = $ (x) + 
n 



o 



L. (x 2 -!, + R 

dx n 



(3.6) 



6/n 



where R = 0 
n 



i/n' 



, and the components are assumed to have 



densities. Here y^ = E 



(X - E[X]) 



and the term 



3 Y 1 



is the conventional dimensionless skewness measure for a dis- 
tribution, being zero for symmetric distributions (Normal), 
and being +2 for the Exponential. Additional terms involving 
kurtosis (4th moments) improve the approximation, but it is 
possible that Edgeworth numerical values can be "infeasible": 
the approximation can actually decrease with x in certain 
ranges. Nevertheless the Edgeworth series has been usefully applied, 
even to unequal component situations, for estimating the loss 



of capacity of an electric utility* see Levy & Kahn [1981]. 

A useful alternative is the method of large deviations , 
Feller [1966], Daniels [1954], and others. The ingenious 
idea is to tilt (or sculpture) the df. components so as to make 
a Normal approximation more effective at predicting the prob- 
ability that S n > x for large x. For equal components with 
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d.f. F(x) look at the tilted probability measure (assumed to 
exist for s > 0, which sometimes restricts the theoretical ap- 
plicability) : 



V (dx) 



e S X F (dx) 
e <Ms) ' 



(3.7) 



/v sX 

'Ms) = £n F (s ) = £n E[e ] being a cumulant generating 
function for F, or X. Manipulations show that 



OO oo 

P{S n > z } = / F n *(dz) = e n ^ (s) l e~ sx v n * (dx) , (3 ' 8) 

z z 

and the idea is to approximate V n * by a Normal centered at z, 
a feat that can be accomplished by choice of s. It turns out 
that it is necessary to solve (sometimes numerically) for s(z) 
the equation 

z = ntj; ' (s ) (3.9) 

in order that the mean of the approximating Normal be at z; 
the variance is n^"(s). Finally, 



P(S n > z} ^ exp n{^[s(z)] - s(z)^'[s(z)] + ~ s 2 (z)^"[s(z)]} * 



1 

/ 2 ? 



oo -—v 

/ e dv 
s ( z) /n tjj" Ls (z) J 



(3.10) 



The above technique can also be applied to a compound 
Poisson model (n is replaced conditionally by N(t), the 
counting process of a Poisson process, and the condition then 
removed) . Such models are frequently employed in inventory 
studies; apparently the large deviations approximation has not 
been applied in that area. 
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(4) Bernoulli Trials and Poisson Process: Rare Events 



Bernoulli Trials are a special case of the Repeated Trials 
model, with events occurring ("Success") or not permitted to 
occur ("Failure") at specific integer time points, often equally 
spaced. In practice the fixed intervals between trials may 
be largely arbitrary, and it is attractive to think of events 

occurring at any (real-valued) time; from this comes the Poisson 
process. One approach to the P.P. properties is to consider a 
B.T. process to operate over time t with unit time steps, and 
then refine the time steps (e.g. let t = 1 day and starting with 
possible demands at 15-minute intervals, then down to 7.5 minutes, 
then to 3.75 etc . ) to create a sequence of B.T. models. The 
limit of the sequence after ultimate refinement describes the 
P.P. 

Specifically, let T(k) be the generic time between successes 

in the (kth) B.T. model with time steps 1/2 (k = 0,1,2, ...). 

This means that the actual number of steps in time t for B.T. 

model k is 2 t; correspondingly, let the probability of success 

k 

per step be P/2 . By conditioning on the first step's outcome 
this means that 



E [ T (k ) ] = l/2 k + (1 - p/2 k ) • E[T (k) ] 



so E[T(k)] = — for every model, as should be true. 
P 

as k -> 00 so time steps become arbitrarily small , 

.k 



P{T(k) > t} = (1 - p/2 k ) fc * 2 -*■ e tp 



(4.1) 
Furthermore , 

(4.2) 



inter-event times become exponentially distributed in the (P.P.) 
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limit. Furthermore the number of P.P. events ("Successes") 
in time t have the Poisson distribution. 

The P.P. is usefully invoked for many modelling purposes. 
Situation . Consider a sequence of days on which demands for 
computer service (time) are made, and focus on the occurrence 
patterns of runs (uninterrupted sequences) of high-deman d days. 
Question: what is the distribution of times between successive 

runs, and what is the distribution of the number of such runs 
in a fixed time t? It will turn out that if either the run 
lengths are long, or if the probability of a high-demand day 
is small, that runs tend to occur as a Poisson process if the 
time scale is appropriate. 

Model . Begin by modelling individual high demand day occurrences 
as successes in B.T. Let x^ (k) represent the time until the 
1st occurrence of a run of length k, and, measured from the end 
of such a run, let (k) , x^ (k) , ... x^ (k) , . . be the time until 

the 2nd, 3rd, ... ith, such run is realized. By the B.T. as- 
sumption (x^(k), i = 1,2, ...} is an IID sequence of rv. 

Then we can represent x (k) by conditioning on the events that 



may occur in the first 


k 


trials 


: 








k 




with 


prob. 


k 

P 




1 


+ x’ (k) 


with 


prob. 


q 


x (k) = - 


j 


+ x' (k) 


with 


prob . 


1 




k 


+ x' (k) 


with 


prob . 


k-1 

P 
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where t' ( k) is an independent replica of any x (k) : the idea is 

that the process starts over once a failure occurs to spoil a run. 
Alternatively, 



x (k) 



k 



R (k) + x ' (k) 



k 

with prob. p 

k 

with prob. 1-p 



where 



P{R(k) 



= j> = 




-1 

k~ ' 



j = 1,2, 



k. 



(4.a) 



(4.5) 



a truncated geometric distribution. From these come the 
generating function of x (k) , and in principle its distribution: 



E 



x (k) 



1 



k k 



z 



qz 



£ 

1- (pz) 
1-pz 




(4.6) 



Differentiation gives the mean 



% V 

P K (1-P) 

the approximation holding if either p -*■ 0 or k In either 

case the run is a rare event. 

While explicit inversion of the expression for E[z T ^] 
is possible by use of partial fractions, the result is quite 
complicated. On the other hand, look for the distribution of 



E[x (k) ] = k + 



1-P 



kp 

1-P } 



x* (k) = x (k) /E[x (k) ] (4.8) 

when E[x (k) ] becomes large. The expectation 
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(4.9) 



E[e - SX * (k) ] = 



-si (k) /E[ t ( k) ] 



can be obtained from the generating function (by putting 
z = exp[-sp q]); next let either p->-0 (rare individual events) 
or k-*-°° (long runs) to find that this transform converges to 
(1 + s) ■*" . Then by the unicity theorem for transforms (Feller [1966]) 
the normalized rv t* ( k) is approximately unit exponentially 
distributed, i.e. 



P{ t (k) < t E[x (k) ]} 'v 1 - e 



(4.10) 



and furthermore the distribution of the number of k-runs in time t, 
N^(t), is approximately Poisson. Deviation from the Poisson 
(indicated by over-variance) may signify that the underlying 
demand generating process is inhomogeneous or cluster-P r c> ne i- n 
time, and that extra facilities are required to reduce backlogs. 
Examination of runs is one way to check the validity of the 
basic modeling assumption of Bernoulli trials. 

Similar limiting arguments simplify other situations in- 
volving rare events that are generated by even more complicated 
processes. See work on first-passage times for combinations of 
random loads by Gaver , Jacobs and Latonche [1981], 
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(5) Markov Models: General Comments 



The basic theory of Markov chains and processes, both in 
discrete and continuous time, is well introduced in standard 
texts such as Feller [1966], Chung [1967], Karlin & Taylor [1975] 
Kleinrock [ 1976] , and needs no systematic coverage, only review 
and illustration. By way of review, recollect the ideas of 
various possible state space definitions: integers, integer 
and real numbers ( "ages" ) , real numbers (e.g. virtual waiting times 
in queues) ; times (index sets) either discrete and 
equally-spaced or imbedded or continuous time; Markov property 
defined by conditional probabilities ("The future is independent 
of the past, given the present"). Carry on to matrix represen- 
tation of the state probabilities after t (0,1,2,...) time 
steps, forward and backward Chapman-Kolmogorov equations, 
generalize to discrete state Markov chain in continuous time 
with exponential sojourns in states, state classification 
emphasizing irreducible chains and transient chains (with at 
least one absorbing barrier) , recurrent events and first-passage 
times and absorption probabilities, generating functions and 
other transforms. 

Simple Markovian assumptions , i.e. that a scalar state 
rv X(t), where t is time or space, is Markov, introduce de- 
pendence in a plausible and tractable manner. Usually it is 
necessary to assume, for example, that the one-step transition 
probabilities (discrete state, discrete time): 

p ± j = P{X(t) = j |X(t-l) = i), (5.D 
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are time-homogeneous in order to obtain explicit neat solutions. 
Analogous assumptions must be made about discrete-state Markov 
processes in continuous time, wherein A^ is the rate of de- 
parture from state i (exponentially distributed sojourn time 
parameter) , and p^ is the corresponding probability of move 
from i to destination j. Of course a known deterministic 
time dependence, involving daily or weekly cycles, and trends 
can be dealt with by numerically multiplying the transition 
probability matrices. 

More irregular changes in process behaviour can be repre- 
sented as the effect of randomly changing external events, or 
random environments for short. In such models the actual 
primary process transition parameters (e.g. p^ j , or A^) 
change in time under the influence of such environmental 
factors as seismic vibration, temperature and humidity, ocean 
sea state, wind speed or other meteorological effects, or 
variations in personnel effectiveness and propensity for errors. 

Random environment models conveniently postulate that environ- 
mental changes induce simple discrete-state Markovian behavior 
on the basic or primary process parameters; of most interest are para- 
meter changes that occur more slowly than do state changes in the 
basic process. 

Markov modelling of real situations usually involves simplifi- 
cations at certain crucial states. Even then, the answers to 
interesting questions may require extensive computing or simu- 
lation. Astute choices of sub-models or component models, 
e.g. the use of "phase-type" distributions for representing 
arrival and departure processes in queues can be of help, as 
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can the recognition (or plausible imposition) and exploitation 
of special structure; see Neuts [1981], 
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(6) Some Markov Process Problems and Models 



Here are some illustrative situations and corresponding 
Markov chain models. 

Situation (Queueing in discrete time). A servicing facility, e.g. 
a computer system or a programming (or other) consultant, or 
a communication channel, experiences single customer arrivals 
in a random fashion; arrivals enter at the discrete times 
0,1, 2, 3, ... only, and service completions occur only at such 

times. Discuss the nature of the delays and backlogs that 
occur. 

Model 1 . Let the probability of a single arrival at time (epoch) 
t be a^ (t) , where i refers to the number present at that 
epoch. Each arrival must wait at least one time period before 
discharge, even if it immediately enters service upon arrival. 

Let d^ (t) be the probability that an arrival that has been 
in service at t actually departs at t+1. Now let X(t) denote 
the number of arrivals in the system who have not yet completed 
service at time t. Model (X(t)} as a Markov chain with the 
following one-step transition probabilities: 

p i i + 1 (t) = P{X(t+l) = i+1 1 X (t) = i} = [l-d i (t) ]a i (t) 

, i > 1 

p. . n (t) = P{X(t+l) = i-1 | X ( t ) = i} = d . (t) [1-a. (t) ] 
i,i-l i 1 

( 6 . 1 ) 

Po,i (t) “ a o (t) 
p o,o (t) ' 1 - a o (t) 
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Pii(t) 



1 - {[l-d i (t) ]a ± (t) + d i (t)[l-a i (t) ]} 



(t) = 0 otherwise; 

If the number in the system is <_ I, so the state space is finite, 
and P T T (t) = 0 and P T , (t) = d (t) then the probability dis- 
tribution of X(t) for any t can be obtained by numerically 
multiplying the one-step transition matrices, o(t) , with elements 
aiven above : 

Pii (t) = P{X(t) - j | X ( 0 ) = i } = element 

in ith row, jth column of 

t 

P(t) = n p(t') (6.2) 

t ' =0 

This can be done especially easily in APL if the process is 
time-homogeneous, i.e. a^ (t) = a^ , d^ (t) = d^ independent of 
elapsed time. Explicit analytical solutions can rarely be 
found for non-time-homogeneous cases, let alone for time homo- 
geneous cases. If they were available, the solutions would 
generally be very complicated and difficult to interpret. 



Model 1-A . Specialize the above to let eu (t) = a>0 and 

d^ (t) = d>0. If the maximum number in the system is I, there 

is a stationary solution; put s = , a = 1-a, d = 1-d: 

ad 



(d-a) d 
dd-aas 1 



TT . 

J 



(d-a) s ^ 
dd-aas 1 



(6.3) 



IT 



I 



(d-a) as 1 
dd-aas 1 
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